Binomial probability
Peter Kamerman and Tory Madden
22 Aug 2018
Question
For each stimulus intenisty, what is the probability that the set of 12 repeated stimuli include zero?
To assess the probability, we modelled the data using the binomial test with a 50% probability of ‘success’ (positive rating arbitrarily chosen as success) at each stimulus intensity for each participant. Significance was assessed at the \(\alpha\) = 0.05 level, using a two-tailed probability.
SPARS A
Import and inspect data
# Import
data <- read_rds('data/SPARS_A.rds')
# Inspect
glimpse(data)## Observations: 1,823
## Variables: 19
## $ PID <chr> "ID01", "ID01", "ID01", "ID01", "ID01", "ID0...
## $ block <chr> "A", "A", "A", "A", "A", "A", "A", "A", "A",...
## $ block_order <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ trial_number <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 1...
## $ intensity <dbl> 3.75, 1.50, 3.25, 1.50, 3.00, 2.75, 1.00, 2....
## $ intensity_char <chr> "3.75", "1.50", "3.25", "1.50", "3.00", "2.7...
## $ rating <dbl> -10, -40, -10, -25, -20, -25, -40, 2, -40, -...
## $ rating_positive <dbl> 40, 10, 40, 25, 30, 25, 10, 52, 10, 40, 54, ...
## $ EDA <dbl> 18315.239, 13904.177, 11543.449, 20542.834, ...
## $ age <dbl> 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, ...
## $ sex <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,...
## $ panas_positive <dbl> 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, ...
## $ panas_negative <dbl> 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, ...
## $ dass42_depression <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ dass42_anxiety <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ dass42_stress <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ pcs_magnification <dbl> 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,...
## $ pcs_rumination <dbl> 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, ...
## $ pcs_helplessness <dbl> 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, ...Binomial test
# Select columns
data %<>%
select(PID, intensity, rating)
# Nest data by PID and stimulus intensity
data_nest <- data %>%
group_by(PID, intensity) %>%
nest()
# Check number of repeats at each stimulus intensity
#data %>%
# group_by(PID, intensity) %>%
# summarise(repeats = n())
# Generate data
data_nest %<>%
# Add probability of success column
mutate(prob = 0.5) %>%
# Extract rating data from dataframe
mutate(data_vec = map(.x = data,
~ .$rating)) %>%
# Recode rating data as categories according to sign
mutate(data_cat = map(.x = data_vec,
~ ifelse(.x < 0,
yes = 'negative',
no = 'positive'))) %>%
# Count the number of positive and negative ratings
## positive numbers arbitrarily listed first == 'success'
mutate(success_count = map(.x = data_cat,
~ c(length(.x[.x == 'positive']),
length(.x[.x == 'negative'])))) %>%
# Conduct binomial test (two-sided)
mutate(binomial_test = map2(.x = success_count,
.y = prob,
~ binom.test(x = .x,
p = .y,
alternative = 'two.sided'))) %>%
# Extract p-value from binomial_test
mutate(binomial_p.value = map(.x = binomial_test,
~ .x$p.value %>%
round(., 3))) %>%
# Categorise p-value using a p < 0.05 threshold
## Significant: distribution deviates significantly
## from the theoretical distribution
## No correction for multiple comparisons
## (too conservative for explorartory analysis)
mutate(significant_p.value = map(.x = binomial_p.value,
~ ifelse(.x < 0.05,
yes = 'yes',
no = 'no')))Plot
For each paticipant, we plotted raw SPARS ratings at each stimulus intensity and colour-coded the data according to whether the p-value returned by the binomial test was significant (distribution of data points deviates significantly from the theoretical expected distribution).
data_nest %>%
# Select data columns
select(PID, intensity, significant_p.value) %>%
# Unnest data
unnest() %>%
# Join with original data
right_join(data) %>%
# Reclass intensity as an ordered factor
mutate(intensity = factor(intensity,
ordered = TRUE)) %>%
# Plot
ggplot(data = .) +
aes(x = intensity,
y = rating,
fill = significant_p.value) +
geom_hline(yintercept = 0) +
geom_point(colour = '#000000',
shape = 21,
size = 3,
alpha = 0.7) +
labs(title = "SPARS A: Binomial test of positive/negative rating distribution",
subtitle = "Probability of 'success' = 0.5* | alpha = 0.05 | two-tailed p-value",
caption = "*'success' arbitrarily chosen as positive SPARS ratings",
x = 'Stimulus intensity (J)',
y = 'SPARS rating') +
scale_x_discrete(breaks = seq(1, 4, by = 0.5),
labels = sprintf('%.2f', seq(1, 4, by = 0.5))) +
scale_y_continuous(limits = c(-50, 50)) +
scale_fill_viridis_d(name = 'Deviates significantly from expected distribution: ',
option = 'C') +
scale_colour_viridis_d(name = 'Deviates significantly from expected distribution: ',
option = 'C') +
facet_wrap(~ PID,
ncol = 4) +
theme(legend.position = 'top',
panel.grid = element_blank(),
panel.spacing = unit(0.1, 'lines'),
strip.text = element_text(margin = margin(t = 0.1,
b = 0.1,
r = 1,
l = 1,
'lines')),
axis.text.x = element_text(angle = -90,
hjust = 0,
vjust = 0.5))
SPARS B
Import and inspect data
# Import
data <- read_rds('data/SPARS_B.rds')
# Inspect
glimpse(data)## Observations: 2,149
## Variables: 10
## $ PID <chr> "FST1601a", "FST1601b", "FST1602", "FST1603", "...
## $ block <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ trial_number <int> 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3,...
## $ intensity <dbl> 3.75, 2.75, 4.25, 3.00, 3.50, 2.50, 2.25, 4.00,...
## $ intensity_rank <int> 7, 3, 8, 3, 5, 4, 1, 8, 3, 3, 8, 3, 6, 4, 4, 8,...
## $ scale <chr> "SPARS", "SPARS", "NRS", "SRS", "SRS", "SPARS",...
## $ scale_combined <chr> "SPARS", "SPARS", "COMBINED", "COMBINED", "COMB...
## $ rating <dbl> 3, 1, 9, -78, -32, 17, 27, 4, 5, 7, 0, -69, 9, ...
## $ rating_eq <dbl> 3.0, 1.0, 4.5, -39.0, -16.0, 17.0, 13.5, 4.0, 5...
## $ rating_z <dbl> 0.5043367, 0.4190808, -0.3963035, -0.6689292, 0...Binomial test
# Select data
data %<>%
# Extract trials rated using the SPARS
filter(scale == 'SPARS') %>%
# Remove ratings of 0
# Ratings of 0 are uninformative
filter(rating != 0) %>%
# Select columns
select(PID, intensity, rating)
# Check number of repeats at each stimulus intensity
data %>%
group_by(PID, intensity) %>%
summarise(repeats = n())## # A tibble: 63 x 3
## # Groups: PID [?]
## PID intensity repeats
## <chr> <dbl> <int>
## 1 FST1601a 2.25 12
## 2 FST1601a 2.5 12
## 3 FST1601a 2.75 12
## 4 FST1601a 3 12
## 5 FST1601a 3.25 8
## 6 FST1601a 3.5 11
## 7 FST1601a 3.75 9
## 8 FST1601a 4 12
## 9 FST1601a 4.25 12
## 10 FST1601b 2.25 11
## # ... with 53 more rows# Nest data by PID and stimulus intensity
data_nest <- data %>%
group_by(PID, intensity) %>%
nest()
# Generate data
data_nest %<>%
# Add probability of success column
mutate(prob = 0.5) %>%
# Extract rating data from dataframe
mutate(data_vec = map(.x = data,
~ .$rating)) %>%
# Recode rating data as categories according to sign
mutate(data_cat = map(.x = data_vec,
~ ifelse(.x < 0,
yes = 'negative',
no = 'positive'))) %>%
# Count the number of positive and negative ratings
## positive numbers arbitrarily listed first == 'success'
mutate(success_count = map(.x = data_cat,
~ c(length(.x[.x == 'positive']),
length(.x[.x == 'negative'])))) %>%
# Conduct binomial test (two-sided)
mutate(binomial_test = map2(.x = success_count,
.y = prob,
~ binom.test(x = .x,
p = .y,
alternative = 'two.sided'))) %>%
# Extract p-value from binomial_test
mutate(binomial_p.value = map(.x = binomial_test,
~ .x$p.value %>%
round(., 3))) %>%
# Categorise p-value using a p < 0.05 threshold
## Significant: distribution deviates significantly
## from the theoretical distribution
## No correction for multiple comparisons
## (too conservative for explorartory analysis)
mutate(significant_p.value = map(.x = binomial_p.value,
~ ifelse(.x < 0.05,
yes = 'yes',
no = 'no')))Plot
For each paticipant, we plotted raw SPARS ratings at each stimulus intensity and colour-coded the data according to whether the p-value returned by the binomial test was significant (distribution of data points deviates significantly from the theoretical expected distribution).
data_nest %>%
# Select data columns
select(PID, intensity, significant_p.value) %>%
# Unnest data
unnest() %>%
# Join with original data
right_join(data) %>%
# Reclass intensity as an ordered factor
mutate(intensity = factor(intensity,
ordered = TRUE)) %>%
# Plot
ggplot(data = .) +
aes(x = intensity,
y = rating,
fill = significant_p.value) +
geom_hline(yintercept = 0) +
geom_point(colour = '#000000',
shape = 21,
size = 3,
alpha = 0.7) +
labs(title = "SPARS B: Binomial test of positive/negative rating distribution",
subtitle = "Probability of 'success' = 0.5* | alpha = 0.05 | two-tailed p-value",
caption = "*'success' arbitrarily chosen as positive SPARS ratings",
x = 'Stimulus intensity (J)',
y = 'SPARS rating') +
scale_x_discrete(breaks = seq(1, 4.5, by = 0.25),
labels = sprintf('%.2f', seq(1, 4.5, by = 0.25))) +
scale_y_continuous(limits = c(-50, 50)) +
scale_fill_viridis_d(name = 'Deviates significantly from expected distribution: ',
option = 'C') +
scale_colour_viridis_d(name = 'Deviates significantly from expected distribution: ',
option = 'C') +
facet_wrap(~ PID,
ncol = 4,
scales = 'free_x') +
theme(legend.position = 'top',
panel.grid = element_blank(),
panel.spacing.x = unit(0.1, 'lines'),
panel.spacing.y = unit(0.5, 'lines'),
strip.text = element_text(margin = margin(t = 0.1,
b = 0.1,
r = 1,
l = 1,
'lines')),
axis.text.x = element_text(angle = -90,
hjust = 0,
vjust = 0.5))
NRS (zero: 0)
Import and inspect data
# Import
data <- read_rds('data/SPARS_B.rds')
# Inspect
glimpse(data)## Observations: 2,149
## Variables: 10
## $ PID <chr> "FST1601a", "FST1601b", "FST1602", "FST1603", "...
## $ block <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ trial_number <int> 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3,...
## $ intensity <dbl> 3.75, 2.75, 4.25, 3.00, 3.50, 2.50, 2.25, 4.00,...
## $ intensity_rank <int> 7, 3, 8, 3, 5, 4, 1, 8, 3, 3, 8, 3, 6, 4, 4, 8,...
## $ scale <chr> "SPARS", "SPARS", "NRS", "SRS", "SRS", "SPARS",...
## $ scale_combined <chr> "SPARS", "SPARS", "COMBINED", "COMBINED", "COMB...
## $ rating <dbl> 3, 1, 9, -78, -32, 17, 27, 4, 5, 7, 0, -69, 9, ...
## $ rating_eq <dbl> 3.0, 1.0, 4.5, -39.0, -16.0, 17.0, 13.5, 4.0, 5...
## $ rating_z <dbl> 0.5043367, 0.4190808, -0.3963035, -0.6689292, 0...Binomial test
# Select data
data %<>%
# Extract trials rated using the SPARS
filter(scale == 'NRS') %>%
# Select columns
select(PID, intensity, rating)
# Check number of repeats at each stimulus intensity
#data %>%
# group_by(PID, intensity) %>%
# summarise(repeats = n())
# Nest data by PID and stimulus intensity
data_nest <- data %>%
group_by(PID, intensity) %>%
nest()
# Generate data
data_nest %<>%
# Add probability of success column
mutate(prob = 0.5) %>%
# Extract rating data from dataframe
mutate(data_vec = map(.x = data,
~ .$rating)) %>%
# Recode rating data as categories according to whether
# the value is greater than 0 (minimum rating on NRS)
mutate(data_cat = map(.x = data_vec,
~ ifelse(.x == 0,
yes = 'negative',
no = 'positive'))) %>%
# Count the number of positive and negative ratings
## positive numbers arbitrarily listed first == 'success'
mutate(success_count = map(.x = data_cat,
~ c(length(.x[.x == 'positive']),
length(.x[.x == 'negative'])))) %>%
# Conduct binomial test (two-sided)
mutate(binomial_test = map2(.x = success_count,
.y = prob,
~ binom.test(x = .x,
p = .y,
alternative = 'greater'))) %>%
# Extract p-value from binomial_test
mutate(binomial_p.value = map(.x = binomial_test,
~ .x$p.value %>%
round(., 3))) %>%
# Categorise p-value using a p < 0.05 threshold
## Significant: distribution deviates significantly
## from the theoretical distribution
## No correction for multiple comparisons
## (too conservative for explorartory analysis)
mutate(significant_p.value = map(.x = binomial_p.value,
~ ifelse(.x < 0.05,
yes = 'yes',
no = 'no')))Plot
For each paticipant, we plotted raw SPARS ratings at each stimulus intensity and colour-coded the data according to whether the p-value returned by the binomial test was significant (distribution of data points deviates significantly from the theoretical expected distribution).
data_nest %>%
# Select data columns
select(PID, intensity, significant_p.value) %>%
# Unnest data
unnest() %>%
# Join with original data
right_join(data) %>%
# Reclass intensity as an ordered factor
mutate(intensity = factor(intensity,
ordered = TRUE)) %>%
# Plot
ggplot(data = .) +
aes(x = intensity,
y = rating,
fill = significant_p.value) +
geom_hline(yintercept = 0) +
geom_point(colour = '#000000',
shape = 21,
size = 3,
alpha = 0.7) +
labs(title = "NRS (zero: 0): Binomial test of positive/negative rating distribution",
subtitle = "Probability of 'success' = 0.5* | alpha = 0.05 | one-tailed (greater) p-value",
caption = "*'success' chosen as NRS rating values greater than 0",
x = 'Stimulus intensity (J)',
y = 'NRS rating') +
scale_x_discrete(breaks = seq(1, 4.5, by = 0.25),
labels = sprintf('%.2f', seq(1, 4.5, by = 0.25))) +
scale_y_continuous(limits = c(0, 100)) +
scale_fill_viridis_d(name = 'Deviates significantly from expected distribution: ',
option = 'C') +
scale_colour_viridis_d(name = 'Deviates significantly from expected distribution: ',
option = 'C') +
facet_wrap(~ PID,
ncol = 4,
scales = 'free_x') +
theme(legend.position = 'top',
panel.grid = element_blank(),
panel.spacing.x = unit(0.1, 'lines'),
panel.spacing.y = unit(0.5, 'lines'),
strip.text = element_text(margin = margin(t = 0.1,
b = 0.1,
r = 1,
l = 1,
'lines')),
axis.text.x = element_text(angle = -90,
hjust = 0,
vjust = 0.5))
NRS (zero: 0 to 15)
Import and inspect data
# Import
data <- read_rds('data/SPARS_B.rds')
# Inspect
glimpse(data)## Observations: 2,149
## Variables: 10
## $ PID <chr> "FST1601a", "FST1601b", "FST1602", "FST1603", "...
## $ block <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ trial_number <int> 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3,...
## $ intensity <dbl> 3.75, 2.75, 4.25, 3.00, 3.50, 2.50, 2.25, 4.00,...
## $ intensity_rank <int> 7, 3, 8, 3, 5, 4, 1, 8, 3, 3, 8, 3, 6, 4, 4, 8,...
## $ scale <chr> "SPARS", "SPARS", "NRS", "SRS", "SRS", "SPARS",...
## $ scale_combined <chr> "SPARS", "SPARS", "COMBINED", "COMBINED", "COMB...
## $ rating <dbl> 3, 1, 9, -78, -32, 17, 27, 4, 5, 7, 0, -69, 9, ...
## $ rating_eq <dbl> 3.0, 1.0, 4.5, -39.0, -16.0, 17.0, 13.5, 4.0, 5...
## $ rating_z <dbl> 0.5043367, 0.4190808, -0.3963035, -0.6689292, 0...Binomial test
# Select data
data %<>%
# Extract trials rated using the SPARS
filter(scale == 'NRS') %>%
# Select columns
select(PID, intensity, rating)
# Check number of repeats at each stimulus intensity
#data %>%
# group_by(PID, intensity) %>%
# summarise(repeats = n())
# Nest data by PID and stimulus intensity
data_nest <- data %>%
group_by(PID, intensity) %>%
nest()
# Generate data
data_nest %<>%
# Add probability of success column
mutate(prob = 0.5) %>%
# Extract rating data from dataframe
mutate(data_vec = map(.x = data,
~ .$rating)) %>%
# Recode rating data as categories according to whether
# the value is greater than 0 (minimum rating on NRS)
mutate(data_cat = map(.x = data_vec,
~ ifelse(.x <= 15,
yes = 'negative',
no = 'positive'))) %>%
# Count the number of positive and negative ratings
## positive numbers arbitrarily listed first == 'success'
mutate(success_count = map(.x = data_cat,
~ c(length(.x[.x == 'positive']),
length(.x[.x == 'negative'])))) %>%
# Conduct binomial test (two-sided)
mutate(binomial_test = map2(.x = success_count,
.y = prob,
~ binom.test(x = .x,
p = .y,
alternative = 'greater'))) %>%
# Extract p-value from binomial_test
mutate(binomial_p.value = map(.x = binomial_test,
~ .x$p.value %>%
round(., 3))) %>%
# Categorise p-value using a p < 0.05 threshold
## Significant: distribution deviates significantly
## from the theoretical distribution
## No correction for multiple comparisons
## (too conservative for explorartory analysis)
mutate(significant_p.value = map(.x = binomial_p.value,
~ ifelse(.x < 0.05,
yes = 'yes',
no = 'no')))Plot
For each paticipant, we plotted raw SPARS ratings at each stimulus intensity and colour-coded the data according to whether the p-value returned by the binomial test was significant (distribution of data points deviates significantly from the theoretical expected distribution).
data_nest %>%
# Select data columns
select(PID, intensity, significant_p.value) %>%
# Unnest data
unnest() %>%
# Join with original data
right_join(data) %>%
# Reclass intensity as an ordered factor
mutate(intensity = factor(intensity,
ordered = TRUE)) %>%
# Plot
ggplot(data = .) +
aes(x = intensity,
y = rating,
fill = significant_p.value) +
geom_hline(yintercept = 0) +
geom_point(colour = '#000000',
shape = 21,
size = 3,
alpha = 0.7) +
labs(title = "NRS (zero: 0 to 15): Binomial test of positive/negative rating distribution",
subtitle = "Probability of 'success' = 0.5* | alpha = 0.05 | one-tailed (greater) p-value",
caption = "*'success' chosen as NRS rating values greater than 15",
x = 'Stimulus intensity (J)',
y = 'NRS rating') +
scale_x_discrete(breaks = seq(1, 4.5, by = 0.25),
labels = sprintf('%.2f', seq(1, 4.5, by = 0.25))) +
scale_y_continuous(limits = c(0, 100)) +
scale_fill_viridis_d(name = 'Deviates significantly from expected distribution: ',
option = 'C') +
scale_colour_viridis_d(name = 'Deviates significantly from expected distribution: ',
option = 'C') +
facet_wrap(~ PID,
ncol = 4,
scales = 'free_x') +
theme(legend.position = 'top',
panel.grid = element_blank(),
panel.spacing.x = unit(0.1, 'lines'),
panel.spacing.y = unit(0.5, 'lines'),
strip.text = element_text(margin = margin(t = 0.1,
b = 0.1,
r = 1,
l = 1,
'lines')),
axis.text.x = element_text(angle = -90,
hjust = 0,
vjust = 0.5))
Session information
sessionInfo()## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] bindrcpp_0.2.2 magrittr_1.5 forcats_0.3.0 stringr_1.3.1
## [5] dplyr_0.7.6 purrr_0.2.5 readr_1.1.1 tidyr_0.8.1
## [9] tibble_1.4.2 ggplot2_3.0.0 tidyverse_1.2.1
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.18 cellranger_1.1.0 pillar_1.3.0
## [4] compiler_3.5.1 plyr_1.8.4 bindr_0.1.1
## [7] tools_3.5.1 digest_0.6.15 viridisLite_0.3.0
## [10] lubridate_1.7.4 jsonlite_1.5 evaluate_0.11
## [13] nlme_3.1-137 gtable_0.2.0 lattice_0.20-35
## [16] pkgconfig_2.0.1 rlang_0.2.1 cli_1.0.0
## [19] rstudioapi_0.7 yaml_2.2.0 haven_1.1.2
## [22] withr_2.1.2 xml2_1.2.0 httr_1.3.1
## [25] knitr_1.20 hms_0.4.2 rprojroot_1.3-2
## [28] grid_3.5.1 tidyselect_0.2.4 glue_1.3.0
## [31] R6_2.2.2 fansi_0.2.3 readxl_1.1.0
## [34] rmarkdown_1.10 modelr_0.1.2 backports_1.1.2
## [37] scales_0.5.0.9000 htmltools_0.3.6 rvest_0.3.2
## [40] assertthat_0.2.0 colorspace_1.3-2 labeling_0.3
## [43] utf8_1.1.4 stringi_1.2.4 lazyeval_0.2.1
## [46] munsell_0.5.0 broom_0.5.0 crayon_1.3.4